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Abstract 

We compute the condensate in QCD with two flavors of dynamical fermions using numerical 
simulation. The simulations use overlap fermions, and the condensate is extracted by fitting the 
distribution of low lying eigenvalues of the Dirac operator in sectors of fixed topological charge to 
the predictions of Random Matrix Theory. 
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I. INTRODUCTION 



The quark condensate £ = (0|gg|0) is the order parameter associated with the sponta- 
neous breaking of chiral symmetry in QCD. Along with the pseudoscalar decay constant, 
it is one of the two fundamental parameters of the lowest order chiral Lagrangian which 
is the low energy effective field theory for QCD. As such, its value is an interesting physi- 
cal quantity whose determination presents a challenge to lattice QCD methodology. While 
there have been many calculations of £ in quenched QCD, results for QCD with dynamical 
fermions are relatively sparse (for a summary of recent results, see Ref. This paper is 
a computation of £ for QCD with two flavors of dynamical overlap fermionsj^, jj. Overlap 
fermions implement chiral symmetry exactly at nonzero lattice spacing via the Ginsparg- 
Wilsonj^] relation. This paper is an extension of our earlier work published in Ref. jlj], which 
was primarily about algorithms, and Ref. 6|, where we use the same methodology to extract 
E in N f = 1 QCD. 

Rather than measure (qq) directly, we will determine the particular combination of the 
coefficients of the low energy effective field theory, £ = f 2 B, in the usual parameterization 



C 2 = ^Tr(d„Ud^) + B^Tv[M(U + U% 



(1) 



One expects that the quantity (qq) (as computed, for example, in a lattice simulation at 
some quark mass m q and simulation volume V) is a function of £, /, m q , and V. For the 
remainder of this paper, we will refer to the quantity £ as the condensate, and this is the 
quantity we will compute on the lattice. 

This is done using the low-lying eigenvalues of the QCD Dirac operator in a finite vol- 
ume, whose distribution can be predicted by random matrix theory (RMT) . This 



hypothesis has 



ulations 



3een checked extensively by lattice calculations, mainly in quenched sim- 



13, 



3, 11, 12, Q. 



i4i, H, Q, but also in dynamical ones using staggered quarks 



181 ] . Our analysis is based on the distribution of the k-th eigenvalue from RMT as 



presented in Refs. [18J, [19( • The prediction is for the distribution of the k-th. eigenvalue of 
the Dirac operator, in each topological sector, which is a function of the dimensionless 
quantity ( = X^V, where E is the chiral condensate and V is the volume of the box. These 
distributions are universal and depend only on the number of flavors and the topological 
charge. They depend parametrically on the dimensionless quantity m q 'EV. By comparing 
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the distribution of the eigenvalues with the RMT prediction one can thus measure the chiral 
condensate S. This method gives the zero quark mass, infinite volume condensate directly. 

The validity of the approach can be verified by comparing the shape of the distribution 
for the various modes and topological sectors. A too small volume causes deviations in the 
shape, particularly for the higher modes. Two recent large scale studies using the overlap 
operator on quenched configurations f° un d that with a length larger than 1.2 fm and 

1.5 fm respectively, the RMT predictions match the result of the simulation. Our dynamical 
lattices have a spatial extent of about 1.5 fm and we have a smaller volume from our earlier 
work. As we will see, random matrix theory describes our data quite well. 

The RMT predictions are made with the assumption that the volume is infinite. In 
the epsilon-regime of chiral perturbation theory, finite volume modifies the formula for the 
condensate by multiplication by a shape factor, a 3 S — > pa 3 S, where 



and c(k/l) depends on the geometry [2 lj. (It is 0.1405 for hypercubes.) We do not know 
p since we have not measured f n , but combining our lattice spacing and lattice size with 
f n = 93 MeV, we expect p ~ 1.42 for our simulations. In what follows, we will refer to the 
quantity we extract from RMT fits as S rather than the more proper label of pS. Ignoring 
p introduces a systematic uncertainty into our result. 

The lattice calculation has four parts: The first is the algorithm. Our calculations are 
performed with overlap fermions, which possess exact chiral symmetry at nonzero lattice 
spacing. A variation on the standard Hybrid Monte Carlo allows us to perform simulations 



in sectors of fixed topology 22j, |23|, |24| . We review the algorithm in Sec. Ill Al 

The extraction of £ from a fit to the eigenvalues of the Dirac operator is described in 

Sec. una 

Next, we need a lattice spacing to convert the dimensionless lattice- regulated condensate 
to a dimensionful number. We obtain the lattice spacing through the Sommer parameter 
from the static quark potential [25j. This is described in Sec. Ill CI 



Finally, we need a matching factor, to convert the lattice-regulated condensate to its MS 



value. We do this using the Regularization Independent scheme [2 61] . This is described in 
Sec. IHD1 

To anticipate our results, which are summarized in Sec. 1111} the larger simulation vol- 



ume considerably improves the quality of the RMT fits. The value of the lattice-regulated 
condensate we obtain here is quite consistent with our earlier result. We find, however, that 
the nonperturbative matching factor we determine is quite different from the perturbative 
matching factor used in Ref. [5j. Our final answer is X 1 / 3 = 282(10) MeV; the quoted error 
is statistical (from the simulation); in addition, we estimate that this number exceeds the 
true result by a systematic factor p 1 / 3 = 1.13. 

II. METHODOLOGY 

A. Lattice action and simulation algorithm 

The massless overlap {2], Q] operator is 

D = D ov (m = 0) = Rq [1 + l5 e(h{-Ro))) (3) 



where e(h) = h/yh 2 is the sign function of the Hermitian kernel operator h = 75G? which is 
taken at negative mass Rq. The squared Hermitian overlap operator H 2 = (75D) 2 = D^D 
commutes with 75 and therefore can have eigenvectors with definite chirality. The modes at 
zero and 4Rq aside (which are associated with zero modes of D), the spectrum is doubled 
with a positive and a negative chirality eigenvector for eigenvalue |A| 2 . 

The RMT analysis requires data restricted to particular topological sectors. It is very 
convenient to generate these data sets directly, rather than letting the topology vary in the 



simulation and filtering it into different topological sectors. A simple observation 22j, |23|, \2A \ 
allows us to do this: the fermion determinant for each flavor is equal to the determinant of 
H 2 in one chiral sector times a correction factor for the modes at m and 2i? : 

det£> = (m/2R ) lQl detH* (4) 



where H 2 pp i s the squared Dirac operator in the chiral sector without zero modes. Since H, 



is a positive operator, its determinant can be replaced by a pseudofermion estimator [241] 



2 

opp 



involving a single chiral pseudofermion for each quark flavor. Then an ensemble can be 



generated using Hybrid Monte Carlo 27]. 



At topological boundaries, the spectrum of the Dirac operator is discontinuous, and so 
is the fermionic contribution to the action. In the algorithm of Ref. [28], the molecular 
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dynamics trajectory either "reflects" from the boundary, with no topological change, or 
"refracts" and changes its topology. We are going to use the algorithm for simulations on 
sectors of fixed topology, so we simply forbid refractions. We pick chiral sources either in the 
opposite chirality (if Q 7^ 0) or randomly select a source chirality (if Q — 0). For other parts 
of the calculation, where we need data sets which are not restricted to particular topological 
sectors, we use our implementation of the algorithm of Ref. 



We now mention specific features of the simulation. Our particular implementation of 
the Hybrid Monte Carlo algorithm has been previously discussed in Refs. fl, We 



use the Liischer-Weisz gauge action 



3l| with the tadpole improved coefficients of Ref. [32 ]. 
Instead of determining the fourth root of the plaquette expectation value u = (([/ p ;)/3) 1//4 
self-consistently, we set it to 0.86 for all our runs as we did in our previous publications. 
We are using a planar kernel Dirac operator d with nearest and next-to-nearest ( "\/2" ) 
interactions. The choice of coefficients, clover term, and value of Ro are those of Refs. {5. 29. 
3o| . Our kernel operator d is constructed from gauge links to which two levels of isotropic 
stout blocking [33J] have been applied. The blocking parameter p is set to 0.15. The sign 
function is computed using the Zolotarev approximation with an exact treatment of the 
low-lying eigenmodes |A) of h(—Ro)- We removed 16 kernel eigenmodes, which typically 
extended in magnitude up to about 0.25 and allowed us to restrict the range of the Zolotarev 
approximation from 0.9 of the maximum value to 2.7 (the upper limit for eigenvalues of our 
kernel action). 

We simulate on 10 4 lattices at one value of the gauge coupling f3 = 7.2 (which we chose 
to be roughly at the N t = 6 phase transition), with three values of the bare sea quark mass 
am q = 0.015, 0.03 and 0.05. 

At each mass value and for \Q\ = 0, 1 we ran several independent data streams: three 
Q = streams for each mass, three \Q\ = 1 streams at am q = 0.05, and two \Q\ = 1 
streams at am q = 0.03 and 0.015. Each stream ran for about 150 trajectories. We used 
one extra level of Hasenbusch preconditioning 34] and typically broke the time intervals into 
eight steps with the lightest mass, four to six steps with the heavier mass, and 12 steps of 
Sexton- Weingarten integration for the gauge fields. Our acceptance rated averaged about 70 
per cent. The whole simulation consumed about 120 processor-months on our array of 3 Ghz 
P4E processors. We dropped the first 50 trajectories from each stream for thermalization 
(and checked that our results were independent of this selection). We recorded eigenvalues 
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after every fifth HMC trajectory. 

By running in the opposite chirality sector, we nearly eliminated critical slowing down: 
the elapsed time for a trajectory was nearly independent of Q and was only marginally 
slower at am q = 0.015 than at am q = 0.05. This would not have been the case had we used 
a conventional (non-chiral) HMC algorithm, because of the near-zero mode. 

This calculation was about a factor 2.5 times slower than our Nf = 1 simulation described 
in Ref. in terms of the cost of the application of an overlap operator to a trial vector. 
(We did that calculation after we collected the data for this project.) This is due to the 
lower level of stout smearing (two steps with p — 0.15 versus three steps for Nf = 1.) 



B. Data fitting 

We computed the lowest four eigenvalues |A| 2 of the squared Dirac operator H 2 . From 
those we get the eigenvalues of the overlap operator D. These lie on a circle, which we 
project onto the imaginary axis via a Mobius transform 

* = 1-A/W (5) 

The distribution of these eigenvalues is predicted by RMT in terms of one parameter E. 
Because we are at finite volume, however, this prediction is valid only for the lowest modes. 
From which level on the deviation between the observed data and the RMT curves becomes 
significant is a priori not clear. Basically, RMT makes two predictions: The positions of the 
peaks are equally spaced and it also predicts the shape of the distribution of each individual 
level. 

These two predictions should be met by the data in order for the extraction of £ to be 
meaningful. The quality of the fit can be measured by the confidence level given by the 
Kolmogorov-Smirnov test {35I . It compares the integrated distributions (the cumulants) of 
the measured data C(x) and the theoretical prediction P(x). The cumulant of the measured 
data is C(x) = n(x)/N where n(x) is the number of data points with a value smaller than 
x and N the total number of data points. The theoretical prediction for this quantity 
can be computed by integrating the distribution: P(x) = f_ f(y)dy. This eliminates the 
bias introduced by binning the data and comparing it to the theoretical distribution. The 
quantity of interest is the largest deviation of P and C: D = max x . \P(x) —C(x)\. From this 
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the confidence level is given by 



Qks ((ViV + 0.12 + 0. 11 /y/N)D 



) 



(6) 



with the function Qks given by 



CO 



Q^(x) = 2^(-r 1 exp(-2jV). 



(7) 



Let us now start by discussing the am q = 0.05 data because of its larger statistics. In 
Fig. [TJ the solid line shows the measured cumulant for the lowest three modes in v = 
and \u\ = 1. This is to be compared to the integrated RMT prediction, which depends 
on the parameter S. We fit S by maximizing the product over the confidence levels for 
each cumulant included in the fit. It is a priori not clear for how many modes the RMT 
description is valid. Due to the finite volume, the distributions are expected to deviate for 
the higher levels. We therefore try several combinations of the number of modes included 
in the fits which we now discuss. The results of these fits can be found in Table [H 

The result of a combined fit to the lowest mode in each topological sector is shown in 
Fig. [TJ To visualize the range of the uncertainty in a 3 E = 0.0105(4), we show the theoretical 
curves for the two extremal values of the one sigma range. The errors on the fit parameter 
are determined by the bootstrap procedure. (See Fig. [2] (top) for a comparison of the binned 
data with the RMT distribution.) We find very good agreement for the shape and position 
of the lowest three modes in v = 0. In v — 1 however, only the lowest mode matches. For 
all these modes, the confidence levels are above 50%. See Table |T] for the confidence levels 
of the fits to the individual distributions. 

To check for consistency between the two topological sectors, we show the results from 
separate fits of the RMT prediction to the distribution in each topological sector. The 
general statements from the combined fit persist: The lowest three modes in v = and the 
lowest one in \u\ = 1 can be described by RMT. The extracted values of a 3 S agree within 
uncertainties (the v — ensemble is completely independent from the \v\ — 1 ensemble). 
We get a 3 £ = 0.0098(5) from v = and a 3 S = 0.0111(6) from v = 1. 

So far we have focused on the am q = 0.05 data set, because of its higher statistics. 
The streams at am q = 0.03 and 0.015, however, give us the possibility to check, whether 
the quark mass dependence of the RMT prediction is correct. From a combined fit to the 
distribution of the lowest eigenvalue in each topological sector (the fit we also chose to use 
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FIG. 1: The measured cumulants for the lowest three modes in v = (left) and |^| = 1 (right). 
The smooth dashed curve represents the RMT prediction with the value of £ from the combined 
fit to the lowest mode in each sector. The dotted curves indicate the la range of S determined by 
the bootstrap method. 

for our final value for am q = 0.05), we extract a 3 S = 0.0112(7) and a 3 S = 0.0130(7) for 
am q = 0.03 and 0.015 respectively. Fig. [2] shows the RMT predictions from these values of 
S along with the measured distributions. 

The values of a 3 S for the three sea quark masses seem to be only barely compatible, 
however, one has to compare them with the difference in the lattice spacing taken into 
account. Using the results from Sec. Ill CI and Sec. Ill Dl (to convert our numbers to MS 
regularization scheme) we plot function of the bare sea quark mass and obtain a 

reasonable agreement within uncertainties, see Fig. [3j Taking the average of the three values 



The measured confidence levels are much better than what we found on smaller lattices 
which makes us confident that the deviations which we observe are only due to the finite 
volume. In Ref. we simulated on 8 4 lattices with a same lattice spacing. The volume 
there was clearly too small. To demonstrate the improvement due to the larger volume, 



we find r 3 E = 0.368(41). 
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FIG. 2: RMT fit to eigenvalue distributions for am q = 0.05, 0.03 and 0.015. The value of a 3 S is 
determined from a combined fit to the lowest mode in each topological sector. 
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m a 3 S \v\ N 


CL level 1 CL level 2 CL level 3 


26 

0.015 0.0130(7) 

1 32 




0.82 
0.42 


0.86 0.02 
0.15 0.71 


33 

0.03 0.0112(7) 

1 28 




0.12 
0.20 


0.57 0.04 
0.13 0.02 


55 

0.05 0.0105(4) 

1 47 




0.55 
0.64 


0.77 0.94 
0.01 0.00 


55 

0.05 0.0098(5) 

1 47 




0.52 
0.67 


0.79 0.93 
0.02 


55 

0.05 0.0111(6) 

1 47 




0.07 
0.99 


0.03 0.001 
0.08 0.0005 


55 

0.05 0.0106(4) 

1 47 




0.52 


0.80 0.90 


0.67 


0.02 0.00 



TABLE I: a 3 S from a fit to the lowest eigenvalues on the 10 4 ensembles. The confidence levels 
for the individual distributions are given in the last three columns; they are boxed if the level is 
included in the fit. 




0.4 



0.00 0.02 0.04 



0.06 



FIG. 3: The dimensionless quantity TqS in MS scheme at /x = 2 GeV, as a function of the bare 
quark mass. The horizontal lines represent the averaged value r|]S = 0.368(41) and its uncertainty. 
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m a 3 £ \u\ N 


CL level 1 CL level 2 CL level 3 


112 

0.03 0.0124(5) 

1 58 


6 • 10~ 5 
9 • 10~ 5 


0.06 0.05 
0.002 1 • 10" 5 


96 

0.05 0.0132(5) 

1 75 




0.01 
0.04 


0.004 0.03 
0.003 



TABLE II: To demonstrate the effect of the finite volume, we also re-analyzed our old 8 4 data. 
The notation is the same as in Tab. HI This table corresponds to the second and third entry in 
Tab. U for the 10 4 lattices. The confidence levels are considerably lower. 

we have re-analyzed the data with the techniques described above. The results are given 
in Table [Til and should be compared to the second and third entry in Table [TJ because in 
both cases we fitted to the lowest mode in each topological sector. We clearly observe 
an improvement in the confidence level when we go to the larger volume, i.e. the RMT 
prediction match the measured cumulants much better. 

We recall (see Eq. [2]) that the result £ of finite simulation volume is equal to pSoo, where 
Sqo is the infinite volume condensate. We do not know p since we have not measured /„■, 
but combining our lattice spacing and lattice size with f n = 93 MeV, we expect p ~ 1.42 
for L = 10 and p(L = 10)/p(L = 8) = 0.86. In order to compare the results from the two 
volumes, we have to divide each of them by its respective p. For m = 0.05 we therefore 
obtain for the 10 4 lattices a 3 ^ = 0.0074(3) and for the 8 4 lattices a 3 Eoo = 0.0080(3). The 
error-bars of the two values touch. However, we suspect that the 8 4 data may be too small 
for RMT to provide a reliable determination of the condensate. We are unwilling to attempt 
a combined fit of p£ in terms of f n and X from data obtained on several simulation volumes 
without good RMT fits for each volume. 



C. Lattice spacing 

We determined an overall scale from a fit to the static quark potential. It is extracted 
from the effective masses of Wilson loops after one level of HYP smearing where 
the short- distance effects of the HYP smearing are corrected using a fit to the perturbative 
lattice artifacts. We measured the potential on 76, 82, and 48 8 3 x 12 configurations at 
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m q 


r /a 


a[fm] 


0.015 


3.47(8) 


0.144(3) 


0.03 


3.27(6) 


0.153(3) 


0.05 


3.35(7) 


0.148(3) 



TABLE III: The Sommer parameter ro and the lattice spacing a in fm from ro = 0.5 fm. We show 
results from two fit ranges used to extract the potential. The fit range for the fit to the potential 
was r G [1.4,6.1] with little variation between different choices for this range. 

am q = 0.05, 0.03, and 0.015 (spaced 5 HMC trajectories apart, from several streams per 
mass) . To generate these data sets we used a conventional two- flavor HMC algorithm. Fits 
with the minimum distance t — 4 and 5 are consistent within uncertainties. In Ref. we 
only had 8 4 lattices, which prevented us from going to these separations. The determinations 
of ro/a from those data sets are slightly different from what we have here, and were slightly 
contaminated by excited states. A compilation of fit results is shown in Table IIHI 



D. Conversion to MS regularization 

A renormalization constant is needed to convert the lattice result of the quark condensate 



to its MS value. To get this matching factor, we use the RI' scheme introduced in Ref. 



26] 



Specifically we follow the procedure described in Ref. [36j. The RI' scheme result in the 



chiral limit can be converted to the 2 GeV MS value by using the ratio connecting the two 



schemes. The ratio was computed by continuum perturbation theory to three loops 37], l38 |. 

The simulation for getting the matching factor should not be restricted in topological 
sectors. Also, one needs simulations with a momentum scale short enough to be free from 
nonperturbative effects and yet not too short to be affected by discretization effects. Some- 
what to our surprise, our sample of 8 4 lattices proved suitable to do this. We combined 
multiple streams of data together to reduce simulation-time autocorrelations. Our data sets 
consisted of approximately 50 configurations per mass value at two values of the bare sea 
quark mass, am q =0.03 and 0.05. 

The 8 4 lattice is periodic in space directions and antiperiodic in the time direction. There- 
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FIG. 4: for quark mass am q = 0.03 and 0.05. 

fore the momentum values are 

aP>1 = \~8 ~8 y, ~8 8^ 2kt + 1 H ■ ^> 

We choose the values of such that the momentum values lie as close as possible to 
the diagonal of the Brillouin zone. The maximum value of ap = 2.115 corresponds to 
ky, = (2, 1,0, 1). The propagators are cast from a point source and then projected to the 
desired momentum values. 

The renormalization constant for the scalar density in the RT scheme is given in Fig. HI 
Values of Z^ 1 ' for am q = 0.03 and 0.05 are listed in Table HVl The inverse lattice spacings 
are 1.29 GeV and 1.32 GeV for the am q = 0.03 and 0.05 data sets determined from the 
Sommer parameter. Thus \i — 2 GeV corresponds to a/i = 1.55 or 1.52. The 2 GeV RF 
values are obtained from linear interpolations from the two closest /i values of the data. 

The conversion ratio for the scalar and pseudoscalar densities, in Landau gauge and to 
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TABLE IV: Values of Z$ in the RT scheme for the two quark masses. The inverse lattice spacings 
are 1.29 GeV and 1.32 GeV for the am q = 0.03 and 0.05 data sets respectively from the Sommer 
parameter. Therefore /i = 2 GeV corresponds to a\x = 1.55 or 1.52 accordingly. The 2 GeV RT 
values are obtained from linear interpolations from the two closest fi values of the data. 





am q = 


0.03 am, = 0.05 


0.878 


0.66(8) 


0.77(13) 


1.178 


0.67(6) 


0.82(14) 


1.416 


0.68(5) 


0.70(6) 


1.619 


0.71(4) 


0.75(7) 


1.800 


0.69(3) 


0.77(7) 


1.963 


0.70(3) 


0.70(3) 


2.115 


0.70(3) 


0.71(3) 


fi = 2GeV 


0.70(4) 


0.73(6) 



three loops, is [37 
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152Cs\ /a, 
3 y Utt 

224993C 3 



729 



54 



80C 4 n/ 2960C 5 



s / ■ 



(9) 



27 27 3 9 

where n/ is the number of flavors and ( n is the Riemann zeta function evaluated at n. This 
formula is for zero quark mass. 

To get numerical results of the above ratio, we use the coupling constant from the so- 



called "ay" scheme. As in the appendix of Ref. 
the plaquette to the coupling 



5], from the one- loop expression relating 



In -TrU v 
3 p 



(10) 



where W = 
ay(3.32/a) = 
is calculated 



0.366 and q*a = 3.32 for the tree-level Liischer-Weisz action, we obtain 
0.192 and 0.193 for the am q = 0.03 and 0.05 data sets. Then af^(e" 5/6 3.32/a) 
4l| and run to af^(2 GeV) by using (3 = 29/12vr and (3 1 = 230/48tt 2 for two 
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flavor QCD. We find af s {2 GeV)= 0.210 and 0.213 for the am q = 0.03 and 0.05 data sets 
respectively. Then from Eq.Q with rif = 2, we get Zf s jZf = 1.160 and 1.164 accord- 
ingly. Therefore Zg IS (2 GeV)= 0.81(5) and 0.85(7) for the two data sets. Because these 
two values have rather large uncertainty, we perform the chiral limit by fitting to a constant 
and get Zf^{2 GeV)=0.83(4). 

III. RESULTS 

Combining our results for the lattice regulated condensate, the Sommer parameter, and 
the Z-factor, we find r 3 E(MS,/i = 2 GeV) = 0.328(28), 0.325(31), and 0.451(45) at am q = 
0.05, 0.03, and 0.015 respectively. These values are plotted in Fig. [3j The result from each 
quark mass is the parameter of the Chiral Lagrangian. Averaging them, we find 

r 3 E(MS, /i = 2 GeV) = 0.368(41). (11) 

Taking the real- world value for r$ = 0.5 fm, this is 

E(MS, /2 = 2 GeV) = 0.0225(25)GeV 3 (12) 

or 

(E(MS, /i = 2 GeV)) 1/3 = 282(10)MeV. (13) 

All of these analyses ignore p. If we knew f n , we could divide it out. This would reduce 
E by a factor 1/1.42 if f n = 93 MeV. Since we have no lattice determination of f n from our 
simulation, p represents a systematic excess, which we must quote separately, 31 MeV for 
E 1 / 3 , for example. 

IV. CONCLUSIONS 

A summary of previous calculations of the condensate has recently been given by 
McNeile [l|]. Comparing the summary of quenched determinations there, a three-flavor pre- 
diction by McNeile, our recent N f = 1 result [6] of (S(M,/i = 2 GeV)) 1 / 3 = 0.269(9) GeV 
(which also does not include p) and our Nf = 2 result quoted above, the condensate seems 
to be a quantity which is not very Nf dependent. Our result is also quite consistent with 
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the Gell-Mann-Oakes-Renner relation, the pion mass and (non-lattice) phenomenological 



estimates of the up and down quark masses 



42]. 



This project suffers from a too-small simulation volume, relatively low statistics and we 
just have a single coarse lattice spacing. The small volume is reflected in the limited range 
of eigenmodes which are well-predicted by RMT and by the factor p which is present as 
a systematic in our final answer. (Recall, p approaches unity like 1/L 2 ). In principle, one 
could determine f v by some other method and scale p out of the RMT fit. We also have the 
usual problem of a lattice simulation at one value of the lattice spacing: scaling violations 
are unknown. In hindsight, the three masses we did here are overkill; a single small quark 
mass suffices to give £ from RMT. 

In contrast, the strength of this calculation is the use of a chiral lattice fermion, which 
allows for the preservation of chiral symmetry in the bare action and insures that spontaneous 
chiral symmetry breaking (and explicit chiral symmetry breaking induced by a quark mass) 
occurs at finite lattice spacing exactly as in the continuum. From a simulation point of 
view, critical slowing down turned out to be largely eliminated due to the possibility to 
run in the opposite chirality sector. The definition of the topological charge by the index 
of the Dirac operator allows us to easily monitor the topological charge. This makes it 
possible to produce ensembles at fixed topology which is convenient for the comparison to 
the predictions of RMT, which are at definite topological charge. 
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